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Abstract 

Provided a special function of one variable and some of its derivatives can be accurately computed 
over a finite range, a method is presented to build a series of polynomial approximations of the 
function with a defined relative error over the whole range. This method is easy to implement and 
makes possible fast computation of special functions. 
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I. INTRODUCTION 

It is often necessary to compute with a high precision special functions of one variable 
within a finite range of values. This task can be very difficult and can require a great 
computational time if the function is known, for instance, by an integral representation or 
by a very long expansion. Such functions can be evaluated with a very high precision by 
symbolic manipulation languages, but this is not a very practical method if you need to 
perform calculations in a Fortran code for instance. 

The idea of the method presented here is to compute the function considered and some of 
its derivatives for a special set of points within the range of interest. This can be performed 
by any mean: symbolic manipulation languages or usual computational codes. The relative 
accuracy required for the function determines completely the number of points and their 
positions within the finite range. Once this set of points is calculated, the function at any 
value within the interval can be computed with the required relative accuracy using only the 
information about the function at the point immediately below and the point immediately 
above the value. This is possible by computing a polynomial whose values and values of 
some of its derivatives are equal to the corresponding values for the function to interpolate, 
for the pair of successive points. 

II. INTERPOLATION WITH FIRST DERIVATIVE 

Let us assume that we know exactly a function F and its first derivative F' at two 
points X\ and x 2 . We can easily determine the third degree polynomial P(x) such that 
P( Xl ) = P(x 2 ) = F(x 2 ), P'(zi) = F'(xi), and P'(x 2 ) = F'(x 2 ). The coefficients of 

the interpolating polynomial can be determined by solving a Vandermonde-like system 
but such a system can be quite ill-conditioned. It is preferable to compute directly P(x) 
by a Lagrange-like formula 2]. Actually, the polynomial P(x) which satisfies the conditions 
above is simply given by 
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provided the spline polynomials / and g are characterized by the boundary properties given 
in Table |U The expressions (jAl|) of these spline functions are given in the Appendix. 



TABLE I: Boundary properties of the spline functions / and g for a third degree interpolating 
polynomial. 
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It is possible to estimate the error made by using P(x) instead of F(x) within the interval 
[xi,^]- To simplify calculations, we can perform a translation of the coordinate system in 
order to fix X\ = and F(xi) = 0, and a rotation to get F'(xi) = 1, for instance. If we note 
X2 = h, F(x2) = y and F'ixz) = z, the interpolating polynomial P(x) is given by 

P(x) = x + ^(3y~zh-2h)x 2 + ^ {-2y + zh + h) x 3 . (2) 

With the same conventions, the limited Taylor expansion of the function F around x% = 
is written 

F(x) =x + ^P-x 2 + ^P-x 3 + ^P-x 4 + 0(x 5 ). (3) 
2 o 24 

Computed in x = X2 = h, the expression above and its first derivative give 

, F"(0)^ 2 F W (0) L3 F( 4 )(0), 4 
F(h) =y &h+ ~^h 2 + — ^h 3 + -^^h\ 

F'(h) = z * 1 + F"(0)h + ^-h 2 + ^-p±h 3 , (4) 

if we neglect contributions of higher order terms. We can solve this system to calculate 
F"(0) and F"'(Q) as a function of h, y, z and F^ 4 )(0). We can then replace these two values 
in Eq. (j3J). Using Eq. (J2J), we finally find 

F(x) - P(x) « ^^x 2 (x - hf . (5) 

The function x 2 (x — h) 2 is represented on Fig. Qfor h = 1. Within the interval [0,h], it 
presents only one maximum at x = h/2, and decreases monotonically from this maximum 
toward zero at x = and x = h. It is then possible to evaluate the maximum error within 
the interval [0, h]. Returning to the first notations, we find 

max \F(x) - P(x) | « fa - x 2 f , (6) 
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FIG. 1: Functions x 2 (x — l) 2 and x 3 (x — l) 3 . In both cases, the extremum appearing within the 
interval [0, 1] is located at x = 0.5. 

the maximal error being located near the middle of the interval. 

For a given set of points for which F(x), F'(x) and are known, it is then possible 
to build an interpolating polynomial for each interval and to estimate the error within each 
interval. But it is possible to use Eq. © in a more clever way. Let us assume that you need 
an approximation of a function F within an interval [a, b] with a fixed relative precision e. 
If you can compute F{x), F'(x) and F^'(x) for arbitrary values x within this range, you 
can start from x\ — a to determine a point X2 in such a way that the relative accuracy of 
the interpolating polynomial defined by Eq. is around e within [xi,^]. Then, you can 
calculate a point x% from x 2 i n a similar way, and so on. The general relation is 
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Finally, a point xn+i > b is reached. With the N + 1 triplets (xj, F(xi), F'(xi)), you can 
build, using Eq. ((TJ, a polynomial approximation of F on iV intervals with iV different 
polynomials Pi(x) of the third degree such that 
F(x) - PAx 



F{x) 
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If you want to compute F(x) with x within the range [a, b], you have to localize first the 
interval [xj, Xj + i] which contains x. Then the calculation at x of the third degree interpolating 



polynomial P{(x) within this interval will give the evaluation of F(x) with a relative error 
of e. These two operations can be performed very fast [l|. 



III. INTERPOLATION WITH FIRST AND SECOND DERIVATIVES 



If you can compute higher order derivatives of the function F, you can build better 
polynomial approximations. The fifth degree polynomial P(x) such that P{x\) = F{x\), 
P(x 2 ) = F{x 2 ), P'(x 1 ) = F'(xi), P'{x 2 ) = F'(x 2 ), P"{ Xl ) = F"( Xl ), and P"{x 2 ) = F"(x 2 ) 
is given by 

P(x) = F( Xl ) f ( + F(x 2 ) f 1 ' 2 
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x 2 - Xi 

provided the spline polynomials /, g and k are characterized by the boundary properties 
given in Table |TI] The expressions (|A2J) of these spline functions are given in the Appendix. 



TABLE II: Boundary properties of the spline functions /, g and k for a fifth degree interpolating 
polynomial. 
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Using the same procedure as in the previous section, the error between the function and 
the interpolating polynomial (JHJ) within the interval [0, h] is estimated at 

F(x) - P(x) ^ ^^-x 3 (x - h) 3 . (10) 

The function x 3 (x — h) 3 is represented on Fig. ^for h — 1. Within the interval [0, h], it 
also presents only one extremum at x = h/2, and tends monotonically from this extremum 
toward zero at x = and x = h. With the most general notations, we find 

max \F(x) - P(x)\ « (x, - x 2 f , (11) 

[xi,X2] 40UoU 



the maximal error being located near the middle of the interval. If you need an approxima- 
tion of a function F with a relative precision e over a fixed range, and if you can compute 
F(x), F'(x), F"(x) and F^\x) for arbitrary values x within this range, you can define a 
series of points with the following relation 
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X i+ i = Xi+ JW(x~) ' ^ ' 

in such a way that the fifth degree polynomials built with Eq. for each interval are an 
approximation of F with the relative accuracy e. 

It is possible to define better and better polynomial approximations by using higher order 
derivatives of the function under study. But very good results can already be obtained with 
the use of the first and second derivatives only. 



IV. APPLICATION AND CONCLUDING REMARKS 

These techniques are used here to compute an approximation of the modified Bessel 
function of integer order Kq{x) [jj. For a fixed range, the number of points decreases if 
the second derivative is used to compute the approximation, in supplement of the first 
derivative only. It is also possible to reduce the number of points by smoothing the function 
to compute. For instance, we have 

v ( \ ~ /^exp(-x) 

*o(s)«^-^, (13) 

for large values of x. If we remove the rapidly varying exponential part of K (x) by comput- 
ing exp(x) Kq(x), we can reduce strongly the number of intervals. The gain is even better by 
computing an approximation of y/x exp(x) Kq(x). These results are illustrated in Table ITTT1 

In order to remove divergent or rapidly varying behaviors, it is sometimes interesting to 
multiply the function F to approximate by a function G known with a very weak relative 
error. An approximation of F(x) G(x) is then computed. The relative precision of the 
approximation of F is not spoiled by dividing the interpolating polynomial by the function 
G, since the relative error on a quotient is the sum of the relative errors of the factors. So, 
if the relative precision for G(x) is very good, the relative error on F(x) is controlled by the 
relative error on F(x) G(x). 

The number of points necessary to reach a fixed precision obviously increases with the 
required accuracy. It depends also strongly on the range of values. This is shown in Table llVl 
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TABLE III: Number of points necessary to reach a relative precision of 10~ 10 for the function F(x) 
with x within the interval [2,6] (Kq(x) is a modified Bessel function). 

F{x) With first derivative With first and 

second derivatives 

K (x) 342 41 

exp(x)K (x) 121 21 

y/x exp(x) Kq(x) 68 15 



TABLE IV: Number of points necessary to reach a relative precision e with first and second 
derivatives for the function ^fx exp(x) Kq[x) within two intervals (Kq(x) is a modified Bessel 
function). 
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The method used here to compute an approximation of a function F over a finite range 
with a definite precision is useful mainly in two cases: 

• You need a code to compute the function F in an usual programming language, but 
the computation with a high accuracy of the function and some of its derivatives is 
only possible in a symbolic manipulation language. 

• You can compute the function F and some of its derivatives in an usual programming 
language, but the calculation time is prohibitive. This can be the case if F is known 
by an integral representation or by a very long expansion, for instance. 

In both cases, it is interesting to compute and store the numbers Xi, F(xi), F'(xi), etc. to 
build a polynomial approximation of F. A demo program is available via anonymous FTP 
on: ftp : / / ftp . umh . ac . be/pub/ f tp_pnt/interp/ 
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APPENDIX A: SPLINE FUNCTIONS 

We give here the spline functions to define the two kinds of interpolating polynomials 
considered in this paper. A third degree interpolating polynomial is defined with the two 
polynomial spline functions 

f{x) = 2x 3 - 3x 2 + 1, 

g(x) = x 3 - 2x 2 + x. (Al) 

Their boundary properties are given in Table |U A fifth degree interpolating polynomial is 
defined with the three polynomial spline functions 

f( x ) = -6x 5 + 15x 4 - I0x 3 + 1, 
g(x) = — 3x 5 + 8x 4 — 6x 3 + x, 

k(x) = \ (-x 5 + 3x 4 - 3x 3 + x 2 ) . (A2) 
Their boundary properties are given in Table |H] 
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